
import numpy as np
import pandas as pd
import math as math
import cmath as cmath
from scipy.stats import norm
from scipy import integrate
import seaborn as sns
from scipy import stats
import matplotlib.pyplot as plt
import tensorflow as tf
import time
from xgboost import XGBClassifier
import warnings
from sklearn.exceptions import DataConversionWarning
from pandas.core.common import SettingWithCopyWarning
warnings.filterwarnings(action='ignore', category=DataConversionWarning)
warnings.simplefilter(action="ignore", category=SettingWithCopyWarning)


# characteristic function
def f(phi, kappa, theta, sigma, rho, v0, r, T, s0, status):
    
    a = kappa * theta
    x = math.log(s0)
    
    # remind that lamda is zero
    if status == 1:
        u = 0.5
        b = kappa - rho * sigma
    else:
        u = -0.5
        b = kappa
    
    d = cmath.sqrt((rho * sigma * phi * 1j - b)**2 - sigma**2 * (2 * u * phi * 1j - phi**2))
    g = (b - rho * sigma * phi * 1j + d) / (b - rho * sigma * phi * 1j - d)
    
    C = r * phi * 1j * T + (a / sigma**2)*((b - rho * sigma * phi * 1j + d) * T - 2 * cmath.log((1 - g * cmath.exp(d * T))/(1 - g)))
    D = (b - rho * sigma * phi * 1j + d) / (sigma**2) * ((1 - cmath.exp(d * T)) / (1 - g * cmath.exp(d * T)))
    
    return cmath.exp(C + D * v0 + 1j * phi * x)


# P1 and P2
def p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 , K, status):
    
    
    integrand = lambda phi: (cmath.exp(-1j * phi * cmath.log(K)) * f(phi, kappa, \
                              theta, sigma, rho, v0, r, T, s0, status) / (1j * phi)).real 
    
    return (0.5 + (1 / math.pi) * integrate.quad(integrand, 0, 100)[0]) # approximate indefinite intergral with a definite one

def p1(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K):
    return p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K, 1)

def p2(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K):
    return p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K, 2)


# call price
def call_price(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K):
    
    P1 = p1(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K)
    P2 = p2(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K)
    
    result = (s0 * P1 - K * math.exp(-r * T) * P2)
    
    if result<0 or k-s0>1: # negative result or initial price very very low with respect to k
        result = 0
        
    return result


# generate the stock price simulation
def stock_price_generator (T, n , m, r, S0, k, V0, sigma, theta, kappa, rho, generator, entire_stock=False ):
    dt = T / n
    
    # Brownian motions:
    dw_v = generator.normal(size=(m, n)) * np.sqrt(dt)
    dw_i = generator.normal(size=(m, n)) * np.sqrt(dt)

    dw_s = rho * dw_v + np.sqrt(1.0 - rho ** 2) * dw_i

    # Perform time evolution 
    s = np.empty((m, n + 1)) # initialisation stock prices vector
    s[:, 0] = S0

    v = np.ones(m) * V0

    for t in range(n):
        dv = kappa * (theta - v) * dt + sigma * np.sqrt(v) * dw_v[:, t]
        ds = r * s[:, t] * dt + np.sqrt(v) * s[:, t] * dw_s[:, t]

        v = np.clip(v + dv, a_min=0.0, a_max=None)
        s[:, t + 1] = s[:, t] + ds
        
    if entire_stock==False: # return only the stocks simulation at T
        result = s[:,-1]; 
    else:
        result = s; # return all the stock dynamic
        
    return result


# find the stock price from ST as expected return
def find_expected_payoff(stock_path, k, r, t_max):
    payoff = max(stock_path[-1] - k, 0) # one payoff for each simulation
    c = payoff * np.exp(-r * t_max)     # in case r=0, this step is useless
    
    return c


generator = np.random.default_rng()


# set parameters and constants
r = 0
k = 2
V0 = 0.01
sigma = 0.6
theta = 0.02
kappa = 6.21
T1 = 10
rho = -0.5

S0 = np.arange(0.5,4,0.01)


time_maturity = 3
num_simulations = 100
n = 250 # numbers of division of the time


#df1 = pd.read_csv('datasets/df1.csv')
#df1 = df1.drop('Unnamed: 0', axis=1)

df1 = pd.DataFrame(columns=['stock_zero', 'price'])


# fill the dataset
for s0 in S0:
    s = stock_price_generator (time_maturity, n, num_simulations, r, s0, k, V0, sigma, theta, kappa, rho, generator, entire_stock=True)

    for stock_path in s: # for each simulation
        price = find_expected_payoff(stock_path, k, r, time_maturity)
        
        new_row = {
                      'price':price, 
                      'stock_zero':stock_path[T1],
                  }
        #append row to the dataframe
        df1 = df1.append(new_row, ignore_index=True)


df1.head()


len(df1)

70000


#df1.to_csv('datasets/df1.csv')


# plot one possible stock evolution as example example
stock_evolution = stock_price_generator (time_maturity, n, 1, r, 2, k, V0, sigma, theta, kappa, rho, generator, entire_stock=True)
stock_evolution = stock_evolution[0]
#plt.plot(stock_evolution)
x = np.linspace(0, time_maturity, num=n+1)

plt.plot(x, stock_evolution)
plt.title("Example of stock evolution following the Heston model")
plt.xlabel("time")
plt.ylabel("St")
#plt.savefig('plots/stock_evolution.pdf')
plt.show()


# function which finds the expected payoff 
def find_mc_price(sT, k, r, T):
    final = sT - k
    payoff = [max(x,0) for x in final] # one payoff for each simulation
    payoff = np.asarray(payoff)
    
    expected_payoff = payoff.mean() # mean
    
    c = expected_payoff * np.exp(-r * T)     # in case r=0, this step is useless
    
    return c


num_sim_mc = 1000


df1_mc = pd.DataFrame(columns=['stock_zero', 'mc_price', 'num_price'])


#df1_mc = pd.read_csv('datasets/df1_mc.csv')
#df1_mc = df1_mc.drop('Unnamed: 0', axis=1)


# find the mc price
for s0 in S0:
    sT = stock_price_generator (time_maturity, n, num_sim_mc, r, s0, k, V0, sigma, theta, kappa, rho, generator)

    mc = find_mc_price(sT, k, r, time_maturity)
    
    num_price = call_price(kappa, theta, sigma, rho, V0, r, time_maturity, s0, k)

    new_row = {
                  'mc_price':mc, 
                  'stock_zero':s0,
                  'num_price': num_price,
              }
    #append row to the dataframe
    df1_mc = df1_mc.append(new_row, ignore_index=True)


df1_mc.head()


len(df1_mc)

350


s_entire = stock_price_generator (time_maturity, n, 1000, r, s0, k, V0, sigma, theta, kappa, rho, generator, entire_stock=True)
expected_price = s_entire[:, -1].mean()

# plot plt.figure(figsize=(8,6))
# noinspection PyTypeChecker
ax_lines = plt.axes()

# Make the line plots
t = np.linspace(0, time_maturity, num=n + 1)
ns = 40
for i in range(ns):
    ax_lines.plot(t, s_entire[i, :], lw=1.0)
ax_lines.set(xlabel='Years', ylabel='St',title='Price Simulations')
ax_lines.set_xlim((0, time_maturity))
ax_lines.set_ylim((1, 7))

# Add mean value to line plots
ax_lines.plot([0.0, 1.0], [expected_price, expected_price], lw='2', ls="--", label='mean')
#plt.savefig('plots/mc_simualtion.pdf')
plt.show()


sns.set_style('darkgrid')
sns.distplot(s_entire[:, -1])
p = norm.pdf(x, mu, std)
plt.plot(x, p, 'k')
plt.xlabel('St')
plt.ylabel('Number')
plt.title('Stock Prices distribution')
plt.legend(["pdf from samples", "Fitted normal distribution"])
#plt.savefig('plots/normal_distribution.pdf')
plt.show()


def stock_price_generator_av (T, n ,m, r, S0, k, V0, sigma, theta, kappa, rho, separate=False):
    dt = T / n
    
    m = int(m/2) # with av we can do half of the simulations
    n = int(n)
    
    # Brownian motions:
    dw_v = generator.normal(size=(m, n)) * np.sqrt(dt)
    dw_i = generator.normal(size=(m, n)) * np.sqrt(dt)

    dw_s = rho * dw_v + np.sqrt(1.0 - rho ** 2) * dw_i

    # Perform time evolution 
    s = np.empty((m, n + 1)) # initialisation stock prices vector
    s[:, 0] = S0
    
    s_ant = np.empty((m, n + 1))
    s_ant[:, 0] = S0

    v = np.ones(m) * V0
    v_ant = np.ones(m) * V0

    for t in range(n):
        dv = kappa * (theta - v) * dt + sigma * np.sqrt(v) * dw_v[:, t]
        dv_ant = kappa * (theta - v_ant) * dt + sigma * np.sqrt(v_ant) * dw_v[:, t]
        ds = r * s[:, t] * dt + np.sqrt(v) * s[:, t] * dw_s[:, t]
        ds_ant = r * s_ant[:, t] * dt + np.sqrt(v_ant) * s_ant[:, t] * dw_s[:, t]

        v = np.clip(v + dv, a_min=0.0, a_max=None)
        v_ant = np.clip(v_ant + dv_ant, a_min=0.0, a_max=None)
        
        s[:, t + 1] = s[:, t] + ds
        s_ant[:, t + 1] = s_ant[:, t] + ds_ant
        
    if separate==True:
        return s, s_ant
        
    return np.concatenate((s, s_ant))[:,-1]


mc_price_av = []


# find the mc price
for s0 in S0:
    ST_ant = stock_price_generator_av (time_maturity, n, num_sim_mc, r, s0, k, V0, sigma, theta, kappa, rho, generator)
    mc_av = find_mc_price(ST_ant, k, r, time_maturity)
    
    mc_price_av.append(mc_av)


#df1_mc['mc_price_av'] = mc_price_av


df1_mc.head()


# df1_mc.to_csv('datasets/df1_mc.csv')


df1_mc = pd.read_csv('datasets/df1_mc.csv')
df1_mc = df1_mc.drop('Unnamed: 0', axis=1)
df1_mc.head()


sns.scatterplot(data=df1, x="stock_zero", y="price", alpha =0.3)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='numerical sol')
plt.plot(S0, df1_mc.mc_price, color='black', linewidth=2, label='MC price')
plt.plot(S0, df1_mc.mc_price_av, color='blue', linewidth=2, label='MC price with av')
plt.legend()
plt.xlabel('S0')
plt.ylabel('Option price')
plt.title("Monte Carlo method")
#plt.savefig('plots/mc_normal_av.pdf')
plt.show()


# show lower volatility for av:
s_ant = stock_price_generator_av (time_maturity, n, num_sim_mc, r, 1, k, V0, sigma, theta, kappa, rho, generator)
s = stock_price_generator (time_maturity, n, num_sim_mc, r, 2, k, V0, sigma, theta, kappa, rho, generator)
print(f"Standard deviation of ST for MC: {np.std(s)}")
print(f"Standard deviation of ST for MC with AV: {np.std(s_ant)}")

Standard deviation of ST for MC: 0.24515159047808166
Standard deviation of ST for MC with AV: 0.24376037892958272


start = time.time()
for s0 in S0:
    ST_ant = stock_price_generator (time_maturity, n, num_sim_mc, r, s0, k, V0, sigma, theta, kappa, rho, generator)
    mc_av = find_mc_price(ST_ant, k, r, time_maturity)
end = time.time()
time_cm = end - start
print(f"Time needed for Monte Carlo: {time_cm}")

Time needed for Monte Carlo: 23.60903525352478


start = time.time()
for s0 in S0:
    ST_ant = stock_price_generator_av (time_maturity, n, num_sim_mc, r, s0, k, V0, sigma, theta, kappa, rho, generator)
    mc_av = find_mc_price(ST_ant, k, r, time_maturity)
end = time.time()
time_cmav = end - start
print(f"Time needed for Monte Carlo with AV: {time_cmav}")

Time needed for Monte Carlo with AV: 21.19796919822693


perc_faster = (time_cm - time_cmav) / time_cm *100
perc_faster = round(perc_faster,2)
print(f"Monte Carlo with AV is {perc_faster}% faster tham the usual Monte Carlo")

Monte Carlo with AV is 10.21% faster tham the usual Monte Carlo


num_sim_mc1 = 100
num_sim_mc2 = 5000


mc_price_av1 = []

# find the mc price
for s0 in S0:
    ST_ant = stock_price_generator_av (time_maturity, n, num_sim_mc1, r, s0, k, V0, sigma, theta, kappa, rho, generator)
    mc_av = find_mc_price(ST_ant, k, r, time_maturity)
    
    mc_price_av1.append(mc_av)


mc_price_av2 = []

# find the mc price
for s0 in S0:
    ST_ant = stock_price_generator_av (time_maturity, n, num_sim_mc2, r, s0, k, V0, sigma, theta, kappa, rho, generator)
    mc_av = find_mc_price(ST_ant, k, r, time_maturity)
    
    mc_price_av2.append(mc_av)


df1_mc['mc_price_av1'] = mc_price_av1
df1_mc['mc_price_av2'] = mc_price_av2


plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='numerical sol')
plt.plot(S0, df1_mc. mc_price_av1, color='orange', linewidth=2, label='num_sim = 100', alpha=0.8)
plt.plot(S0, df1_mc. mc_price_av, color='blue', linewidth=2, label='num_sim = 1000')
plt.plot(S0, df1_mc. mc_price_av2, color='magenta', linewidth=2, label='num_sim = 5000')
plt.title("Monte Carlo pricing with AV using different number of simulations")
plt.legend()
plt.xlabel('S0')
plt.ylabel('Option price')
#plt.savefig('plots/mc_different_sim.pdf')
plt.show()


def root_mean_square_error(num_sol, proposed_solution):
    n = len(num_sol) 
    difference = num_sol-proposed_solution
    square_diff = [x ** 2 for x in difference]
    error = np.sqrt((sum(square_diff)/n))
    return error


err_mc = root_mean_square_error(df1_mc.num_price, df1_mc.mc_price)
err_mc_av = root_mean_square_error(df1_mc.num_price, df1_mc.mc_price_av)
err_mc_av1 = root_mean_square_error(df1_mc.num_price, df1_mc.mc_price_av1)
err_mc_av2 = root_mean_square_error(df1_mc.num_price, df1_mc.mc_price_av2)

print(f"RMSQ Monte Carlo (num_simulations = {num_sim_mc}): {err_mc}")
print(f"RMSQ Monte Carlo with AV (num_simulations = {num_sim_mc1}): {err_mc_av1}")
print(f"RMSQ Monte Carlo with AV (num_simulations = {num_sim_mc}): {err_mc_av}")
print(f"RMSQ Monte Carlo with AV (num_simulations = {num_sim_mc2}): {err_mc_av2}")

RMSQ Monte Carlo (num_simulations = 1000): 0.14124623047806095
RMSQ Monte Carlo with AV (num_simulations = 100): 0.15578851031457136
RMSQ Monte Carlo with AV (num_simulations = 1000): 0.1383395065212146
RMSQ Monte Carlo with AV (num_simulations = 5000): 0.1392209677705778


from sklearn.linear_model import LinearRegression, LogisticRegression, Ridge
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import mean_squared_error
from sklearn.metrics import confusion_matrix
from sklearn.tree import DecisionTreeRegressor
from sklearn.metrics import roc_auc_score
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import roc_curve
from sklearn.neighbors import KNeighborsClassifier
from xgboost.sklearn import XGBRegressor
from sklearn import metrics
from sklearn.metrics import classification_report

from keras.metrics import BinaryAccuracy, AUC, Precision, Recall
from keras.models import Sequential
from keras.layers import Dense
from keras.callbacks import History
from keras.optimizers import Adam


df1.head()


plt.hist(df1.price, bins =20)
plt.title("Number of occurances for each price")
#plt.savefig('plots/occurance_price_all.pdf')
plt.show()


sns.boxplot(x=df1.price)
plt.title("Boxplot for price distribution")
#plt.savefig('plots/volatility_prices.pdf')
plt.show()


y = df1.price
X = df1.stock_zero

x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.2)


x_train = x_train.values.reshape(-1, 1)
y_train = y_train.values.reshape(-1, 1)
x_test = x_test.values.reshape(-1, 1)


lin = LinearRegression()
lin.fit(x_train, y_train)
y_pred_lin = lin.predict(x_test)

print("Linear Regression RMSE:", np.sqrt(mean_squared_error(y_test, y_pred_lin)))

Linear Regression RMSE: 0.5615332819145643


price_linear = lin.predict(S0.reshape(-1, 1))


sns.scatterplot(data=df1, x="stock_zero", y="price", alpha =0.3)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='numerical sol')
plt.plot(S0, price_linear, color='black', linewidth=2, label='linear regression')
plt.legend()
#plt.savefig('plots/linear_regression_only.pdf')
plt.show()


df1_class = pd.DataFrame()
df1_class['exercise'] = df1.price!=0
df1_class['stock_zero'] = df1.stock_zero
df1_class.head()


X_class = df1_class['stock_zero']
y_class = df1_class['exercise']

x_train_class, x_test_class, y_train_class, y_test_class = train_test_split(X_class, y_class, test_size=0.2)


x_train_class = x_train_class.values.reshape(-1, 1)
y_train_class = y_train_class.values.reshape(-1, 1)
x_test_class = x_test_class.values.reshape(-1, 1)


logreg = LogisticRegression()
logreg.fit(x_train_class, y_train_class)

LogisticRegression()


execute_log = logreg.predict(S0.reshape(-1, 1))


df1_ml = pd.DataFrame()


price_log_lin = []

for i in range(len(S0)):
    if execute_log[i]==0:
        price_log_lin.append(0)
    else:
        price_log_lin.append(price_linear[i])


sns.scatterplot(data=df1, x="stock_zero", y="price", alpha =0.3)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='numerical sol')
plt.plot(S0, price_linear, color='black', linewidth=2, label='linear regression')
plt.plot(S0, price_log_lin, color='blue', linewidth=2, label='logistic class+linear regr')
plt.legend()
#plt.savefig('plots/linear_regress_class.pdf')
plt.show()


df1_class['price'] = df1.price


df1_class['exercise_predict'] = logreg.predict(df1_class['stock_zero'].values.reshape(-1, 1))


df1_class.head()


# execute linear regression only on nonzero values (only if exercised)
# => liear regression only on cleaned data

df1_class_linear = df1_class.loc[ df1_class['exercise_predict']==True ]
df1_class_linear.head()


y = df1_class_linear.price
X = df1_class_linear.stock_zero

x_train_lin_cleaned, x_test_lin_cleaned, y_train_lin_cleaned, y_test_lin_cleaned = train_test_split(X, y, test_size=0.2)


x_train_lin_cleaned = x_train_lin_cleaned.values.reshape(-1, 1)
y_train_lin_cleaned = y_train_lin_cleaned.values.reshape(-1, 1)
x_test_lin_cleaned = x_test_lin_cleaned.values.reshape(-1, 1)


lin_clean = LinearRegression()
lin_clean.fit(x_train_lin_cleaned, y_train_lin_cleaned)
y_pred_lin_cleaned = lin_clean.predict(x_test_lin_cleaned)

print("Linear Regression RMSE:", np.sqrt(mean_squared_error(y_test_lin_cleaned, y_pred_lin_cleaned)))

Linear Regression RMSE: 0.7056372381920196


price_pred_lin_cleaned = lin_clean.predict(S0.reshape(-1, 1))


df1_ml = pd.DataFrame()


price_log_lin_cleaned = []

for i in range(len(S0)):
    if execute_log[i]==0:
        price_log_lin_cleaned.append(0)
    else:
        price_log_lin_cleaned.append(price_pred_lin_cleaned[i])


plt.figure(figsize=(10,7))
df1_noout = df1[(np.abs(stats.zscore(df1['price'])) < 3)]
a = sns.scatterplot(data=df1_noout, x="stock_zero", y="price", alpha =0.3)
axes = a.axes
axes.set_ylim(-0.25,2.5)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='numerical sol')
plt.plot(S0, df1_mc. mc_price_av, color='y', linewidth=2, label='MC price with av')
plt.plot(S0, price_linear, color='black', linewidth=2, label='linear regression')
plt.plot(S0, price_log_lin, color='blue', linewidth=2, label='logistic class+linear regr')
plt.plot(S0, price_log_lin_cleaned, color='m', linewidth=2, label='logistic class+linear cleaned regr')
plt.legend()
#plt.savefig('plots/linear_regress_class_mc.pdf')
plt.show()


def plot_confusion_matrix (confusion_matrix, title):

    ax = plt.subplot()
    sns.heatmap(confusion_matrix,annot=True,fmt = "d",square = True,ax = ax, linewidths = 1, linecolor = "w", cmap = "Pastel2")
    ax.set_xlabel('True labels')
    ax.set_ylabel('Predicted labels')
    ax.xaxis.set_ticklabels(['Exercised','Not Exercised'])
    ax.yaxis.set_ticklabels(['Exercised','Not Exercised'], va="center")
    plt.title(title)
    fig_directory = 'plots/confusion_matrix_' + title[:5] + '.pdf'
    plt.savefig(fig_directory)
    plt.show()


def plot_roc_curve(X_test, y_test, y_pred, model, title):
    logit_roc_auc = roc_auc_score(y_test, y_pred)
    fpr, tpr, thresholds = roc_curve(y_test, logreg.predict_proba(X_test)[:,1])
    plt.figure()
    plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc)
    plt.plot([0, 1], [0, 1],'r--')
    plt.xlim([0.0, 1.0])
    plt.ylim([0.0, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    titlee = 'Exercise ROC curve ' + title
    plt.title(titlee)
    plt.legend(loc="lower right")
    fig_directory = 'plots/roc_curve_' + title[:5] + '.pdf'
    plt.savefig(fig_directory)
    plt.show()


df1_class = pd.DataFrame()
df1_class['exercise'] = df1.price!=0
df1_class['stock_zero'] = df1.stock_zero
df1_class['binary'] = [1 if x==True else 0 for x in df1_class['exercise']]
df1_class.head()


X = df1_class['stock_zero']
y = df1_class['exercise']

X_train, X_test, y_train, y_test, price_train, price_test = train_test_split(X, y, df1['price'], test_size=0.20)


X_train = X_train.values.reshape(-1, 1)
y_train = y_train.values.reshape(-1, 1)
X_test = X_test.values.reshape(-1, 1)


start_log = time.time()
logreg = LogisticRegression()
logreg.fit(X_train, y_train)
end_log = time.time()
end_log-start_log

0.08766770362854004


y_pred_logistic = logreg.predict(X_test)


confusion_matrix_logistic = confusion_matrix(y_test, y_pred_logistic, labels=[True, False])
plot_confusion_matrix (confusion_matrix_logistic, 'Logistic Regression')


print(classification_report(y_test, y_pred_logistic))

              precision    recall  f1-score   support

       False       0.87      0.86      0.87      6444
        True       0.89      0.89      0.89      7556

    accuracy                           0.88     14000
   macro avg       0.88      0.88      0.88     14000
weighted avg       0.88      0.88      0.88     14000


plot_roc_curve(X_test, y_test, y_pred_logistic, logreg, 'Logistic Reg')


Neighbor_List=[3,5,10,20]
parameters = {'n_neighbors':Neighbor_List}
KNNC = KNeighborsClassifier()
classifier_knn = GridSearchCV(KNNC, parameters, cv=5, verbose=0, scoring ='accuracy')

classifier_knn.fit(X_train, y_train)

GridSearchCV(cv=5, estimator=KNeighborsClassifier(),
             param_grid={'n_neighbors': [3, 5, 10, 20]}, scoring='accuracy')


y_pred_knn = classifier_knn.predict(X_test)


print('Accuracy KNN: {:.2f}%'.format(logreg.score(X_test, y_test)*100))

Accuracy KNN: 87.71%


confusion_matrix_knn = confusion_matrix(y_test, y_pred_knn, labels=[True, False])
plot_confusion_matrix (confusion_matrix_knn, 'KNN')


plot_roc_curve(X_test, y_test, y_pred_knn, classifier_knn, 'KNN')


print(classification_report(y_test, y_pred_knn))

              precision    recall  f1-score   support

       False       0.89      0.84      0.86      6444
        True       0.87      0.91      0.89      7556

    accuracy                           0.88     14000
   macro avg       0.88      0.87      0.87     14000
weighted avg       0.88      0.88      0.88     14000


params = {
    'n_estimators': [100, 200, 400],
    'max_depth': [2, 3, 6],
    'learning_rate': [0.005, 0.01, 0.02],
    'subsample': [0.4, 0.8]
}

classifier_xgboost = GridSearchCV(XGBClassifier(random_state=10), params, scoring ='accuracy')
classifier_xgboost.fit(X_train, y_train)

GridSearchCV(estimator=XGBClassifier(base_score=None, booster=None,
                                     colsample_bylevel=None,
                                     colsample_bynode=None,
                                     colsample_bytree=None, gamma=None,
                                     gpu_id=None, importance_type='gain',
                                     interaction_constraints=None,
                                     learning_rate=None, max_delta_step=None,
                                     max_depth=None, min_child_weight=None,
                                     missing=nan, monotone_constraints=None,
                                     n_estimators=100, n_jobs=None,
                                     num_parallel_tree=None, random_state=10,
                                     reg_alpha=None, reg_lambda=None,
                                     scale_pos_weight=None, subsample=None,
                                     tree_method=None, validate_parameters=None,
                                     verbosity=None),
             param_grid={'learning_rate': [0.005, 0.01, 0.02],
                         'max_depth': [2, 3, 6],
                         'n_estimators': [100, 200, 400],
                         'subsample': [0.4, 0.8]},
             scoring='accuracy')


y_pred_xgboost = classifier_xgboost.predict(X_test)


print(classification_report(y_test, y_pred_xgboost))

              precision    recall  f1-score   support

       False       0.90      0.83      0.86      6444
        True       0.87      0.92      0.89      7556

    accuracy                           0.88     14000
   macro avg       0.88      0.88      0.88     14000
weighted avg       0.88      0.88      0.88     14000


confusion_matrix_knn = confusion_matrix(y_test, y_pred_xgboost, labels=[True, False])
plot_confusion_matrix (confusion_matrix_knn, 'XGBoost')


plot_roc_curve(X_test, y_test, y_pred_xgboost, classifier_xgboost, 'XGBoost')


def print_recall(train_metric, val_metric):
    # summarize history for loss
    plt.plot(train_metric)
    plt.plot(val_metric)
    plt.plot([0, len(train_metric)], [0.5, 0.5], '--', color='black', alpha=0.7)
    plt.title('Model recall progress')
    plt.ylabel('accuracy')
    plt.xlabel('epoch')
    plt.ylim(ymax = 1, ymin = 0)
    plt.legend(['train', 'test'], loc='upper left')
    plt.show()


epoch_num = 30


# 1 hidden layer
start_nn = time.time()
model = Sequential()
model.add(Dense(20, input_dim=1, activation='relu'))
model.add(Dense(10, activation='relu'))
model.add(Dense(1, activation='sigmoid'))

model.compile(optimizer=Adam(),
              loss='binary_crossentropy',
              metrics=[BinaryAccuracy(), AUC(), Precision(), Recall()])


start_nn = time.time()
history = model.fit(X_train, y_train, epochs=epoch_num, validation_data=(X_test, y_test), batch_size=128, callbacks=[History()])
end_nn = time.time()

Epoch 1/30
438/438 [==============================] - 2s 4ms/step - loss: 0.5392 - binary_accuracy: 0.6922 - auc_1: 0.8704 - precision_1: 0.6486 - recall_1: 0.9298 - val_loss: 0.3664 - val_binary_accuracy: 0.8776 - val_auc_1: 0.9465 - val_precision_1: 0.8474 - val_recall_1: 0.9367
Epoch 2/30
438/438 [==============================] - 1s 2ms/step - loss: 0.3083 - binary_accuracy: 0.8776 - auc_1: 0.9464 - precision_1: 0.8765 - recall_1: 0.8983 - val_loss: 0.2984 - val_binary_accuracy: 0.8784 - val_auc_1: 0.9466 - val_precision_1: 0.8702 - val_recall_1: 0.9046
Epoch 3/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2877 - binary_accuracy: 0.8778 - auc_1: 0.9485 - precision_1: 0.8766 - recall_1: 0.8987 - val_loss: 0.2903 - val_binary_accuracy: 0.8801 - val_auc_1: 0.9465 - val_precision_1: 0.8634 - val_recall_1: 0.9180
Epoch 4/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2836 - binary_accuracy: 0.8789 - auc_1: 0.9486 - precision_1: 0.8696 - recall_1: 0.9107 - val_loss: 0.2876 - val_binary_accuracy: 0.8796 - val_auc_1: 0.9466 - val_precision_1: 0.8644 - val_recall_1: 0.9154
Epoch 5/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2822 - binary_accuracy: 0.8789 - auc_1: 0.9486 - precision_1: 0.8665 - recall_1: 0.9153 - val_loss: 0.2879 - val_binary_accuracy: 0.8794 - val_auc_1: 0.9466 - val_precision_1: 0.8558 - val_recall_1: 0.9277
Epoch 6/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2815 - binary_accuracy: 0.8792 - auc_1: 0.9485 - precision_1: 0.8662 - recall_1: 0.9163 - val_loss: 0.2859 - val_binary_accuracy: 0.8805 - val_auc_1: 0.9466 - val_precision_1: 0.8630 - val_recall_1: 0.9196
Epoch 7/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2807 - binary_accuracy: 0.8793 - auc_1: 0.9486 - precision_1: 0.8647 - recall_1: 0.9187 - val_loss: 0.2894 - val_binary_accuracy: 0.8769 - val_auc_1: 0.9466 - val_precision_1: 0.8443 - val_recall_1: 0.9400
Epoch 8/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2804 - binary_accuracy: 0.8796 - auc_1: 0.9486 - precision_1: 0.8638 - recall_1: 0.9205 - val_loss: 0.2854 - val_binary_accuracy: 0.8792 - val_auc_1: 0.9465 - val_precision_1: 0.8555 - val_recall_1: 0.9279
Epoch 9/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2801 - binary_accuracy: 0.8799 - auc_1: 0.9486 - precision_1: 0.8624 - recall_1: 0.9233 - val_loss: 0.2844 - val_binary_accuracy: 0.8796 - val_auc_1: 0.9466 - val_precision_1: 0.8578 - val_recall_1: 0.9252
Epoch 10/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2799 - binary_accuracy: 0.8797 - auc_1: 0.9486 - precision_1: 0.8623 - recall_1: 0.9231 - val_loss: 0.2846 - val_binary_accuracy: 0.8782 - val_auc_1: 0.9465 - val_precision_1: 0.8523 - val_recall_1: 0.9304
Epoch 11/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2798 - binary_accuracy: 0.8796 - auc_1: 0.9486 - precision_1: 0.8601 - recall_1: 0.9261 - val_loss: 0.2838 - val_binary_accuracy: 0.8801 - val_auc_1: 0.9466 - val_precision_1: 0.8635 - val_recall_1: 0.9179
Epoch 12/30
438/438 [==============================] - 1s 3ms/step - loss: 0.2798 - binary_accuracy: 0.8796 - auc_1: 0.9485 - precision_1: 0.8608 - recall_1: 0.9250 - val_loss: 0.2850 - val_binary_accuracy: 0.8779 - val_auc_1: 0.9466 - val_precision_1: 0.8484 - val_recall_1: 0.9359
Epoch 13/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2796 - binary_accuracy: 0.8794 - auc_1: 0.9486 - precision_1: 0.8599 - recall_1: 0.9260 - val_loss: 0.2836 - val_binary_accuracy: 0.8796 - val_auc_1: 0.9465 - val_precision_1: 0.8578 - val_recall_1: 0.9252
Epoch 14/30
438/438 [==============================] - 1s 3ms/step - loss: 0.2794 - binary_accuracy: 0.8795 - auc_1: 0.9486 - precision_1: 0.8598 - recall_1: 0.9264 - val_loss: 0.2840 - val_binary_accuracy: 0.8784 - val_auc_1: 0.9466 - val_precision_1: 0.8517 - val_recall_1: 0.9318
Epoch 15/30
438/438 [==============================] - 1s 3ms/step - loss: 0.2795 - binary_accuracy: 0.8796 - auc_1: 0.9485 - precision_1: 0.8590 - recall_1: 0.9278 - val_loss: 0.2836 - val_binary_accuracy: 0.8791 - val_auc_1: 0.9465 - val_precision_1: 0.8553 - val_recall_1: 0.9279
Epoch 16/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2796 - binary_accuracy: 0.8796 - auc_1: 0.9485 - precision_1: 0.8595 - recall_1: 0.9270 - val_loss: 0.2837 - val_binary_accuracy: 0.8786 - val_auc_1: 0.9466 - val_precision_1: 0.8524 - val_recall_1: 0.9314
Epoch 17/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2795 - binary_accuracy: 0.8791 - auc_1: 0.9486 - precision_1: 0.8592 - recall_1: 0.9263 - val_loss: 0.2835 - val_binary_accuracy: 0.8796 - val_auc_1: 0.9466 - val_precision_1: 0.8576 - val_recall_1: 0.9257
Epoch 18/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2795 - binary_accuracy: 0.8792 - auc_1: 0.9486 - precision_1: 0.8601 - recall_1: 0.9253 - val_loss: 0.2841 - val_binary_accuracy: 0.8779 - val_auc_1: 0.9466 - val_precision_1: 0.8504 - val_recall_1: 0.9327
Epoch 19/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2796 - binary_accuracy: 0.8798 - auc_1: 0.9486 - precision_1: 0.8600 - recall_1: 0.9267 - val_loss: 0.2835 - val_binary_accuracy: 0.8791 - val_auc_1: 0.9466 - val_precision_1: 0.8554 - val_recall_1: 0.9279
Epoch 20/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2796 - binary_accuracy: 0.8791 - auc_1: 0.9485 - precision_1: 0.8578 - recall_1: 0.9286 - val_loss: 0.2836 - val_binary_accuracy: 0.8804 - val_auc_1: 0.9466 - val_precision_1: 0.8621 - val_recall_1: 0.9207
Epoch 21/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2795 - binary_accuracy: 0.8790 - auc_1: 0.9486 - precision_1: 0.8586 - recall_1: 0.9271 - val_loss: 0.2846 - val_binary_accuracy: 0.8781 - val_auc_1: 0.9466 - val_precision_1: 0.8489 - val_recall_1: 0.9357
Epoch 22/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2794 - binary_accuracy: 0.8790 - auc_1: 0.9486 - precision_1: 0.8596 - recall_1: 0.9256 - val_loss: 0.2837 - val_binary_accuracy: 0.8784 - val_auc_1: 0.9466 - val_precision_1: 0.8523 - val_recall_1: 0.9310
Epoch 23/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2796 - binary_accuracy: 0.8792 - auc_1: 0.9485 - precision_1: 0.8594 - recall_1: 0.9263 - val_loss: 0.2852 - val_binary_accuracy: 0.8776 - val_auc_1: 0.9465 - val_precision_1: 0.8468 - val_recall_1: 0.9379
Epoch 24/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2794 - binary_accuracy: 0.8795 - auc_1: 0.9486 - precision_1: 0.8589 - recall_1: 0.9277 - val_loss: 0.2835 - val_binary_accuracy: 0.8791 - val_auc_1: 0.9466 - val_precision_1: 0.8552 - val_recall_1: 0.9281
Epoch 25/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2793 - binary_accuracy: 0.8793 - auc_1: 0.9486 - precision_1: 0.8586 - recall_1: 0.9277 - val_loss: 0.2841 - val_binary_accuracy: 0.8784 - val_auc_1: 0.9466 - val_precision_1: 0.8496 - val_recall_1: 0.9350
Epoch 26/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2794 - binary_accuracy: 0.8794 - auc_1: 0.9486 - precision_1: 0.8591 - recall_1: 0.9270 - val_loss: 0.2872 - val_binary_accuracy: 0.8756 - val_auc_1: 0.9466 - val_precision_1: 0.8397 - val_recall_1: 0.9446
Epoch 27/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2794 - binary_accuracy: 0.8791 - auc_1: 0.9486 - precision_1: 0.8586 - recall_1: 0.9271 - val_loss: 0.2856 - val_binary_accuracy: 0.8769 - val_auc_1: 0.9465 - val_precision_1: 0.8443 - val_recall_1: 0.9400
Epoch 28/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2795 - binary_accuracy: 0.8792 - auc_1: 0.9485 - precision_1: 0.8587 - recall_1: 0.9273 - val_loss: 0.2835 - val_binary_accuracy: 0.8798 - val_auc_1: 0.9465 - val_precision_1: 0.8589 - val_recall_1: 0.9241
Epoch 29/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2796 - binary_accuracy: 0.8784 - auc_1: 0.9485 - precision_1: 0.8582 - recall_1: 0.9263 - val_loss: 0.2845 - val_binary_accuracy: 0.8785 - val_auc_1: 0.9466 - val_precision_1: 0.8696 - val_recall_1: 0.9055
Epoch 30/30
438/438 [==============================] - 1s 2ms/step - loss: 0.2797 - binary_accuracy: 0.8792 - auc_1: 0.9485 - precision_1: 0.8596 - recall_1: 0.9260 - val_loss: 0.2834 - val_binary_accuracy: 0.8794 - val_auc_1: 0.9466 - val_precision_1: 0.8557 - val_recall_1: 0.9279


end_nn - start_nn

32.66271686553955


print_recall( history.history['recall_1'], history.history['val_recall_1'])


positive_x = X_test[(logreg.predict(X_test)).flatten() > 0.5].flatten()
negative_x = X_test[(logreg.predict(X_test)).flatten() < 0.5].flatten()
positive_y = price_test[(logreg.predict(X_test)).flatten() > 0.5]
negative_y = price_test[(logreg.predict(X_test)).flatten() < 0.5]

sns.scatterplot(x=positive_x, y=positive_y, color='b', alpha =0.3, label="executed")
sns.scatterplot(x=negative_x, y=negative_y, color='r', alpha =0.3, label="not-executed")
plt.title("Classification results with Logistic Regression")
plt.xlabel('S0')
plt.ylabel('Option price')
#plt.savefig('plots/scatter_log.PNG')
plt.show()


positive_x = X_test[(model.predict(X_test)).flatten() > 0.5].flatten()
negative_x = X_test[(model.predict(X_test)).flatten() < 0.5].flatten()
positive_y = price_test[(model.predict(X_test)).flatten() > 0.5]
negative_y = price_test[(model.predict(X_test)).flatten() < 0.5]

sns.scatterplot(x=positive_x, y=positive_y, color='b', alpha =0.3, label="executed")
sns.scatterplot(x=negative_x, y=negative_y, color='r', alpha =0.3, label="not-executed")
plt.title("Classification results with NN")
plt.xlabel('S0')
plt.ylabel('Option price')
#plt.savefig('plots/scatter_nn.PNG')
plt.show()


# we get s0 and the numerical solution from the mc dataset
df1_class_s0 = pd.DataFrame()
df1_class_s0 = df1_mc[['stock_zero', 'num_price']]
df1_class_s0['real_exercise'] =  [False if x==0.0 else True for x in df1_class_s0.num_price]
df1_class_s0 = df1_class_s0.drop('num_price', axis=1)


df1_class_s0['log_regre'] = logreg.predict(df1_class_s0['stock_zero'].values.reshape(-1, 1))
df1_class_s0['knn'] = classifier_knn.predict(df1_class_s0['stock_zero'].values.reshape(-1, 1))
df1_class_s0['xg_boost'] = classifier_xgboost.predict(df1_class_s0['stock_zero'].values.reshape(-1, 1))
df1_class_s0.head()


# from False/True to 0/1 
df1_class_s0 = df1_class_s0.astype(int)
df1_class_s0.head()


error_classification = {}
for column in ['log_regre', 'knn', 'xg_boost']:
    error_classification[column] = root_mean_square_error(df1_class_s0['real_exercise'], df1_class_s0[column])


error_classification

{'log_regre': 0.31622776601683794,
 'knn': 0.27255405754769874,
 'xg_boost': 0.2449489742783178}


df1['pred_executed_logistic'] = logreg.predict(df1['stock_zero'].values.reshape(-1, 1))
df1.head()


df1_regr = df1[df1['pred_executed_logistic']==True]
df1_regr = df1_regr.drop('pred_executed_logistic', axis=1)
df1_regr.head()


plt.hist(df1_regr.price, bins =20)
plt.title("Number of occurances for each price after classification")
#plt.savefig('plots/occurance_price_postclass.pdf')
plt.show()


X = df1_regr['stock_zero']
y = df1_regr['price']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20)


X_train = X_train.values.reshape(-1, 1)
y_train = y_train.values.reshape(-1, 1)
X_test = X_test.values.reshape(-1, 1)


lin = LinearRegression()
lin.fit(X_train, y_train)
y_pred_lin = lin.predict(X_test)


parameters = {'alpha': np.arange(0.001, 0.01, 0.1)}
ridge = Ridge()
ridge_best = GridSearchCV(ridge, parameters, scoring ='neg_mean_squared_error')
ridge_best.fit(X_train, y_train)
y_pred_ridge = ridge_best.predict(X_test)


DTSplitList=[0.001,0.05,0.1, 0.99]
DTLeafList=[0.005, 0.1,0.5]
criterion =['mse']
max_depth = list(range(30,0,-2))
parameters = {'criterion':criterion,'max_depth':max_depth, 'min_samples_leaf': DTLeafList, 'min_samples_split':DTSplitList}
DTR = DecisionTreeRegressor(random_state=10)
dec_tree = GridSearchCV(DTR, parameters, scoring ='neg_mean_squared_error')

dec_tree.fit(X_train, y_train)
y_pred_dec = dec_tree.predict(X_test)


print("DTR best parameters:", dec_tree.best_params_)

DTR best parameters: {'criterion': 'mse', 'max_depth': 30, 'min_samples_leaf': 0.005, 'min_samples_split': 0.1}


RFEstimatorList = [25,50,100,200]
criterion =['mse']
max_depth = list(range(20,0,-2))

parameters = {'criterion':criterion,'max_depth':max_depth,'n_estimators':RFEstimatorList}

RFR = RandomForestRegressor(random_state=10)
RFR_best = GridSearchCV(RFR, parameters, cv=5, verbose=0,scoring ='neg_mean_squared_error')
RFR_best.fit(X_train, y_train)

y_pred_random = RFR_best.predict(X_test)


print("RFR best parameters:", RFR_best.best_params_)

RFR best parameters: {'criterion': 'mse', 'max_depth': 4, 'n_estimators': 200}


params = {
    'n_estimators': [100, 200, 300, 400],
    'max_depth': [2, 3, 6],
    'learning_rate': [0.005, 0.01, 0.02],
    'subsample': [0.4, 0.6, 0.8],
    'objective' :['reg:squarederror']
    
}
xgb_regress = GridSearchCV(XGBRegressor(silent=False,random_state=10), params, scoring ='neg_mean_squared_error')
xgb_regress.fit(X_train, y_train)
y_pred_xgb_regress = xgb_regress.predict(X_test)

print("XGBR best parameters:", xgb_regress.best_params_)


df1_regre = pd.DataFrame(columns=['stock_zero'])


lin_s0 = lin.predict(S0.reshape(-1, 1))
ridge_best_s0 = ridge_best.predict(S0.reshape(-1, 1))
dec_tree_s0 = dec_tree.predict(S0.reshape(-1, 1))
RFR_best_s0 = RFR_best.predict(S0.reshape(-1, 1))
xgb_regress_s0 = xgb_regress.predict(S0.reshape(-1, 1))


df1_regre = pd.read_csv('plots/df1_regre.csv')


prediction_nn_s0 = model_reg.predict(S0.reshape(-1, 1))
df1_regre['nn_regress_s0'] = prediction_nn_s0 
df1_regre.head()


# classification
S0_class = logreg.predict(S0.reshape(-1, 1))


df1_regre['stock_zero'] = S0
df1_regre['linear_best_s0'] = ridge_best_s0
df1_regre['ridge_best_s0'] = ridge_best_s0
df1_regre['dec_tree_s0'] = dec_tree_s0
df1_regre['RFR_best_s0'] = RFR_best_s0
df1_regre['xgb_regress_s0'] = xgb_regress_s0
df1_regre['nn_regress_s0'] = prediction_nn_s0 = model.predict(S0.reshape(-1, 1))


for i in range(len(S0_class)): 
    if S0_class[i]==False:
        df1_regre.iloc[i] = (S0[i], 0, 0, 0, 0, 0, 0, 0)


df1_regre.head()


print("Linear RMSE:", np.sqrt(mean_squared_error(df1_mc.num_price, df1_regre.linear_best_s0)))
print("Ridge RMSE:", np.sqrt(mean_squared_error(df1_mc.num_price, df1_regre.ridge_best_s0)))
print("Decision RMSE:", np.sqrt(mean_squared_error(df1_mc.num_price, df1_regre.dec_tree_s0)))
print("Random RMSE:", np.sqrt(mean_squared_error(df1_mc.num_price, df1_regre.RFR_best_s0)))
print("XGBoost RMSE:", np.sqrt(mean_squared_error(df1_mc.num_price, df1_regre.xgb_regress_s0)))

Linear RMSE: 0.15119602024427806
Ridge RMSE: 0.15119602024427806
Decision RMSE: 0.15398368539366317
Random RMSE: 0.15382917157579842
XGBoost RMSE: 0.1565113853649121


#df1_regre.to_csv('plots/df1_regre.csv')


plt.figure(figsize=(10,7))
#sns.scatterplot(data=df1, x="stock_zero", y="price", alpha =0.3)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='Numerical SOlution')
plt.plot(S0, df1_mc.mc_price_av, color='y', linewidth=2, label='MC with AV')
plt.plot(S0, df1_regre.linear_best_s0, color='orange', linewidth=2, label='Linear Regression')
plt.plot(S0, df1_regre.ridge_best_s0, color='magenta', linewidth=2, label='Ridge Regression')
plt.plot(S0, df1_regre.dec_tree_s0, color='green', linewidth=2, label='Decision Tree')
plt.plot(S0, df1_regre.RFR_best_s0, color='black', linewidth=2, label='Random Forest')
plt.plot(S0, df1_regre.xgb_regress_s0, color='blue', linewidth=2, label='XGBoost')
plt.title("Regression Methods")
plt.legend()
plt.savefig('plots/regression_models.pdf')
plt.show()


df1_ordered = df1.sort_values(by=['stock_zero']).reset_index(drop=True)
df1_ordered.head()


S0_class = logreg.predict(S0.reshape(-1, 1))
X_class = df1_ordered['pred_executed_logistic']


# X:
X_len = len(X_class)
area_size_step = int(X_len/20) # 5% of the total S0 range
area_size_step

3500


# identify step position
where_changes = [i for i in range(1,X_len) if X_class[i]!=X_class[i-1] ]
step_position = int(np.mean(where_changes))
step_position

32525


# where begin regression
begin_regression_pos = step_position - area_size_step
begin_regression_pos

29025


df1_regr_with_area = df1_ordered.loc[begin_regression_pos:]
df1_regr_with_area.head()


X_area = df1_regr_with_area['stock_zero']
y_area = df1_regr_with_area['price']

X_train_area, X_test_area, y_train_area, y_test_area = train_test_split(X_area, y_area, test_size=0.20)


X_train_area = X_train_area.values.reshape(-1, 1)
y_train_area = y_train_area.values.reshape(-1, 1)
X_test_area = X_test_area.values.reshape(-1, 1)


lin_area = LinearRegression()
lin_area.fit(X_train_area, y_train_area)
y_pred_lin_area = lin_area.predict(X_test_area)


# S0:
S0_len = len(S0)
area_size_step_S0 = int(S0_len/10)


where_changes_S0 = [i for i in range(1,len(S0)) if S0_class[i]!=S0_class[i-1] ]
step_position_S0 = int(np.mean(where_changes_S0))
step_position_S0

162


# where begin regression
begin_regression_pos_S0 = step_position_S0 - area_size_step_S0
begin_regression_pos_S0

127


lin_s0_area_regr = lin_area.predict(S0[begin_regression_pos_S0:].reshape(-1, 1))
lin_s0_area_regr = [x if x>0 else 0 for x in lin_s0_area_regr] 
lin_s0_with_area = np.concatenate((np.zeros(begin_regression_pos_S0), lin_s0_area_regr))


plt.figure(figsize=(10,7))
#sns.scatterplot(data=df1, x="stock_zero", y="price", alpha =0.3)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='Numerical SOlution')
plt.plot(S0, df1_mc. mc_price_av, color='y', linewidth=2, label='MC price with av')
plt.plot(S0, df1_regre.linear_best_s0, color='blue', linewidth=2, label='Linear Regression')
plt.plot(S0, lin_s0_with_area, color='black', linewidth=2, label='Linear Regression with Area')
plt.title("Linear regression with area")
plt.legend()
plt.savefig('plots/linear_with_area.pdf')
plt.show()


print("Numerical solution: \n")
print("  Linear RMSE without area:", np.sqrt(mean_squared_error(df1_mc.num_price, df1_regre.linear_best_s0)))
print("  Linear RMSE with area:", np.sqrt(mean_squared_error(df1_mc.num_price, lin_s0_with_area)))
print("\nMonte Carlo: \n")
print("  Linear RMSE without area:", np.sqrt(mean_squared_error(df1_mc. mc_price_av, df1_regre.linear_best_s0)))
print("  Linear RMSE with area:", np.sqrt(mean_squared_error(df1_mc. mc_price_av, lin_s0_with_area)))

Numerical solution: 

  Linear RMSE without area: 0.15119602024427806
  Linear RMSE with area: 0.14213526207369503

Monte Carlo: 

  Linear RMSE without area: 0.0630167939775045
  Linear RMSE with area: 0.0434779198188384


params = {
    'n_estimators': [100, 200, 300],
    'max_depth': [2, 3, 6],
    'learning_rate': [0.01, 0.02],
    'subsample': [0.6, 0.8],
    'objective' :['reg:squarederror']
    
}
tree_area = GridSearchCV(XGBRegressor(silent=False,random_state=10), params, scoring ='neg_mean_squared_error')

tree_area.fit(X_train_area, y_train_area)


tree_area_regr = tree_area.predict(S0[begin_regression_pos_S0:].reshape(-1, 1))
tree_area_regr = [x if x>0 else 0 for x in tree_area_regr] 
tree_with_area = np.concatenate((np.zeros(begin_regression_pos_S0), tree_area_regr))


plt.figure(figsize=(10,7))
#sns.scatterplot(data=df1, x="stock_zero", y="price", alpha =0.3)
plt.plot(S0, df1_mc.num_price, color='r', linewidth=2, label='Numerical SOlution')
plt.plot(S0, df1_mc. mc_price_av, color='y', linewidth=2, label='MC price with av')
plt.plot(S0, df1_regre.xgb_regress_s0, color='blue', linewidth=2, label='XGBoost')
plt.plot(S0, tree_with_area, color='black', linewidth=2, label='XGBoost with area')
plt.title("XGBoost with area")
plt.legend()
plt.savefig('plots/decision_with_area.pdf')
plt.show()


print("Numerical solution: \n")
print("  XGBoost RMSE without area:", np.sqrt(mean_squared_error(df1_mc.num_price,df1_regre.xgb_regress_s0)))
print("  XGBoost RMSE with area:", np.sqrt(mean_squared_error(df1_mc.num_price, tree_with_area)))
print("\nMonte Carlo: \n")
print("  XGBoost RMSE without area:", np.sqrt(mean_squared_error(df1_mc. mc_price_av, df1_regre.xgb_regress_s0)))
print("  XGBoost RMSE with area:", np.sqrt(mean_squared_error(df1_mc. mc_price_av, tree_with_area)))

Numerical solution: 

  XGBoost RMSE without area: 0.1565113853649121
  XGBoost RMSE with area: 0.1512635292718006

Monte Carlo: 

  XGBoost RMSE without area: 0.06450042029852468
  XGBoost RMSE with area: 0.04012366235396293


# model parameters
Sigma = [0.4, 0.8]
Theta = [0.01, 0.05]
Kappa = [5, 6]
Rho = [-0.5, -0.7]
T1 = 10
S0 = np.arange(0.5,4,0.1)
V0 = [0.01, 0.1]


# option parameters
int_rates = [0, 0.05]
k = 2
time_maturity = [1, 3]


num_simulations = 100
n = 100


df2 = pd.DataFrame(columns=['price', 'stock_zero',
                           'sigma',  
                           'maturity', 'initial_vol', 
                           'kappa', 'rho',
                           'theta','interest_rates'])


# fill the dataset
for v0 in V0:
    print(f"\nvo {v0}")
    for sigma in Sigma:
        print(f"sigma {sigma}")
        for theta in Theta:
            print(f"theta {theta}")
            for kappa in Kappa:
                print(f"kappa {kappa}")
                for rho in Rho:
                    print(f"rho {rho}")
                    for s0 in S0:
                        for t in time_maturity:
                            for r in int_rates:
                                s = stock_price_generator (t, n, num_simulations, r, s0, k, v0, sigma, theta, kappa, rho, generator, entire_stock=True)

                                for stock_path in s: # for each simulation
                                    price = find_expected_payoff(stock_path, k, r, t)
                                    
                                    new_row = {
                                        'price':price, 
                                        'stock_zero':stock_path[T1],
                                        'maturity':t, 
                                        'sigma':sigma,
                                        'theta':theta,
                                        'kappa': kappa,
                                        'rho': rho,
                                        'initial_vol': v0, 
                                        'interest_rates':r
                                    }
                                    #append row to the dataframe
                                    df2 = df2.append(new_row, ignore_index=True)

vo 0.01
sigma 0.4
theta 0.01
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7
theta 0.05
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7
sigma 0.8
theta 0.01
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7
theta 0.05
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7

vo 0.1
sigma 0.4
theta 0.01
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7
theta 0.05
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7
sigma 0.8
theta 0.01
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7
theta 0.05
kappa 5
rho -0.5
rho -0.7
kappa 6
rho -0.5
rho -0.7


df2.head()


#df2.to_csv('df2.csv')


df2['exercise'] = df2.price!=0
df2 = df2.drop(['price', 'interest_rates'], axis=1)
df2.head()


df2.groupby('exercise').count().stock_zero

exercise
False    202515
True     245485
Name: stock_zero, dtype: int64


table = pd.crosstab(df2.theta, df2.exercise)
table.div(table.sum(1).astype(float), axis=0).plot(kind='bar', stacked=True)
plt.title('Stacked Bar Chart depending on the theta')
plt.xlabel('theta')
plt.ylabel('Proportion');


table = pd.crosstab(df2.sigma, df2.exercise)
table.div(table.sum(1).astype(float), axis=0).plot(kind='bar', stacked=True)
plt.title('Stacked Bar Chart depending on the sigma')
plt.xlabel('sigma')
plt.ylabel('Proportion');


table = pd.crosstab(df2.initial_vol, df2.exercise)
table.div(table.sum(1).astype(float), axis=0).plot(kind='bar', stacked=True)
plt.title('Stacked Bar Chart depending on the initial vol')
plt.xlabel('initial_vol')
plt.ylabel('Proportion');


table = pd.crosstab(df2.maturity, df2.exercise)
table.div(table.sum(1).astype(float), axis=0).plot(kind='bar', stacked=True)
plt.title('Stacked Bar Chart depending on the maturity')
plt.xlabel('T')
plt.ylabel('Proportion')
plt.savefig('e_vs_ne_strike')


X = df2.drop('exercise', axis = 1)
y = df2['exercise']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20)


X_train = X_train.values.reshape(-1, 1)
X_test = X_test.values.reshape(-1, 1)
y_train = y_train.values.reshape(-1, 1)
y_test = y_test.values.reshape(-1, 1)


logreg2 = LogisticRegression()
logreg2.fit(X_train, y_train)

LogisticRegression()


X.columns

Index(['stock_zero', 'sigma', 'maturity', 'initial_vol', 'kappa', 'rho',
       'theta'],
      dtype='object')


weights = logreg2.coef_
abs_weights = np.abs(weights).flatten()

print(f"Print Feature importance: \n")
for i in range(len(X.columns)):
    print(f"{X.columns[i]}: {abs_weights[i]}")

Print Feature importance: 

stock_zero: 3.3723396053896884
sigma: 0.06440089290982819
maturity: 0.009522686583890211
initial_vol: 0.3605304077959656
kappa: 0.015389850642859334
rho: 0.06372447910342552
theta: 0.05335492


# set parameters and constants
r = 0
K = np.arange(0.5,4,0.05)
V0 = 0.01
sigma = 0.6
theta = 0.02
kappa = 6.21
rho = -0.5

s0 = 2


time_maturity = 3
num_simulations = 100
n = 250 # numbers of division of the time


df3 = pd.DataFrame(columns=['k', 'price'])
num_price = []


# fill the dataset
for k in K:
    s = stock_price_generator (time_maturity, n, num_simulations, r, s0, k, V0, sigma, theta, kappa, rho, generator, entire_stock=True)

    for stock_path in s: # for each simulation
        price = find_expected_payoff(stock_path, k, r, time_maturity)
        num_price.append(call_price(kappa, theta, sigma, rho, V0, r, time_maturity, s0, k))
        
        new_row = {
                      'price':price, 
                      'k':k,
                  }
        #append row to the dataframe
        df3 = df3.append(new_row, ignore_index=True)


df3_class = pd.DataFrame()
df3_class['exercise'] = df3.price!=0
df3_class['k'] = df3.k
df3_class['price'] = df3.price
df3_class.head()


X_class = df3_class['k']
y_class = df3_class['exercise']

x_train_class, x_test_class, y_train_class, y_test_class = train_test_split(X_class, y_class, test_size=0.2)


x_train_class = x_train_class.values.reshape(-1, 1)
y_train_class = y_train_class.values.reshape(-1, 1)
x_test_class = x_test_class.values.reshape(-1, 1)


logreg3 = LogisticRegression()
logreg3.fit(x_train_class, y_train_class)

LogisticRegression()


df3_class['exercise_predict'] = logreg3.predict(df3_class['k'].values.reshape(-1, 1))


df3_class_linear = df3_class.loc[ df3_class['exercise_predict']==True ]
df3_class_linear.head()


y = df3_class_linear.price
X = df3_class_linear.k

x_train_lin_cleaned, x_test_lin_cleaned, y_train_lin_cleaned, y_test_lin_cleaned = train_test_split(X, y, test_size=0.2)


x_train_lin_cleaned = x_train_lin_cleaned.values.reshape(-1, 1)
y_train_lin_cleaned = y_train_lin_cleaned.values.reshape(-1, 1)
x_test_lin_cleaned = x_test_lin_cleaned.values.reshape(-1, 1)


lin_clean3 = LinearRegression()
lin_clean3.fit(x_train_lin_cleaned, y_train_lin_cleaned)

LinearRegression()


price_pred_lin_cleaned3 = lin_clean3.predict(K.reshape(-1, 1))


price_log_lin_cleaned3 = []

for i in range(len(K)):
    if execute_log3[i]==0:
        price_log_lin_cleaned3.append(0)
    else:
        price_log_lin_cleaned3.append(price_pred_lin_cleaned3[i])


num_price = []

for k in K:
    num_price.append(call_price(kappa, theta, sigma, rho, V0, r, time_maturity, s0, k))


plt.plot(K, num_price, color='red', linewidth=2, label='Numerical solution')
plt.plot(K, price_log_lin_cleaned3, color='black', linewidth=2, label='Logistic+Linear') 
plt.legend()
plt.title("Comparison methods when K moves and S0 is fixed")
plt.xlabel("K")
plt.ylabel("price")
plt.show()

	exercise	stock_zero	price	exercise_predict
14188	True	2.124788	0.461163	True
14625	True	2.193899	0.474660	True
14734	False	2.129866	0.000000	True
14825	False	2.133750	0.000000	True
14860	False	2.126172	0.000000	True

	stock_zero	price
14188	2.124788	0.461163
14625	2.193899	0.474660
14706	2.119692	0.256519
14734	2.129866	0.000000
14825	2.133750	0.000000

	stock_zero	price	pred_executed_logistic
29025	1.944814	0.000000	False
29026	1.944823	0.000000	False
29027	1.944867	0.000000	False
29028	1.944921	0.385596	False
29029	1.945016	0.000000	False

Heston Model Project¶

Import packages¶

Explicit solution of the Heston price equation ¶

Stock paths generator ¶

Generate Data ¶

Monte Carlo simulation ¶

Plot MC simulations only to give a graphical idea¶

Monte Carlo simulation with Antithetic variables ¶

Time analysis for Monte Carlo ¶

Changes with different numbers of simulations ¶

Error Analysis MC ¶

Machine Learning models ¶

Motivation in combining regression and classification ¶

Classification methods ¶

Logistic Regression¶

KNN¶

XGBoost Classifier¶

Neural Network¶

Error Classification¶

Regression methods ¶

Linear Regression¶

Ridge Regression¶

Decision Tree¶

Random Forest¶

XGBoost Regression¶

Analyse errors¶

Use different Regression methods with S0¶

ATM area ¶

Add other features ¶

Classification¶

Logistic Regression¶

S0 fixed and K moves ¶

	stock_zero	price
0	0.524303	0.0
1	0.522428	0.0
2	0.499181	0.0
3	0.494353	0.0
4	0.502161	0.0

	stock_zero	mc_price	num_price	mc_price_av
0	0.50	0.0	0.0	0.0
1	0.51	0.0	0.0	0.0
2	0.52	0.0	0.0	0.0
3	0.53	0.0	0.0	0.0
4	0.54	0.0	0.0	0.0

	Unnamed: 0	stock_zero	nn_regress_s0
0	0	0.50	0.420878
1	1	0.51	0.424886
2	2	0.52	0.428893
3	3	0.53	0.432901
4	4	0.54	0.436908

	stock_zero	pred_executed_logistic
0	0.387751	False
1	0.407598	False
2	0.419984	False
3	0.422262	False
4	0.426654	False

	stock_zero	sigma	maturity	initial_vol	kappa	rho	theta	exercise
0	0.515862	0.4	1.0	0.01	5.0	-0.5	0.01	False
1	0.506035	0.4	1.0	0.01	5.0	-0.5	0.01	False
2	0.480199	0.4	1.0	0.01	5.0	-0.5	0.01	False
3	0.503627	0.4	1.0	0.01	5.0	-0.5	0.01	False
4	0.525715	0.4	1.0	0.01	5.0	-0.5	0.01	False

	exercise	k	price
0	True	0.5	1.297998
1	True	0.5	1.637903
2	True	0.5	1.893333
3	True	0.5	0.663243
4	True	0.5	1.571816

	stock_zero	real_exercise	log_regre	knn	xg_boost
0	0	0	0	0	0
1	0	0	0	0	0
2	0	0	0	0	0
3	0	0	0	0	0
4	0	0	0	0	0